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Far-Ultraviolet and Far-Infrared Bivariate Luminosity Function of Galaxies: 
Complex Relation between Stellar and Dust Emission 
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Far-ultraviolet (FUV) and far-infrared (FIR) luminosity functions (LFs) of galaxies show a strong evolution 
from z = to z = 1, but the FIR LF evolves much stronger than the FUV one. The FUV is dominantly radiated 
from newly formed short-lived OB stars, while the FIR is emitted by dust grains heated by the FUV radiation field. 
It is known that dust is always associated with star formation activity. Thus, both FUV and FIR are tightly related 
to the star formation in galaxies, but in a very complicated manner In order to disentangle the relation between 
FUV and FIR emissions, we estimate the UV-IR bivariate LF (BLF) of galaxies with GALEX and AKARl All-Sky 
Survey datasets. Recently we invented a new mathematical method to construct the BLF with given marginals 
and prescribed correlation coefficient. This method makes use of a tool from mathematical statistics, so called 
"copula". The copula enables us to construct a bivariate distribution function from given marginal distributions 
with prescribed correlation and/or dependence structure. With this new formulation and FUV and FIR univariate 
LFs, we analyze various FUV and FIR data with GALEX, Spitzer, and AKARI to estimate the UV-IR BLF. The 
obtained BLFs naturally explain the nonlinear complicated relation between FUV and FIR emission from star- 
forming galaxies. Though the faint-end of the BLF was not well constrained for high-z samples, the estimated 
linear correlation coefficient p was found to be very high, and is remarkably stable with redshifts (from 0.95 at 
z = to 0.85 at z = 1.0). This implies the evolution of the UV-IR BLF is mainly due to the different evolution of 
the univariate LFs, and may not be controlled by the dependence structure. 
Key words: Dust; galaxies: formation; galaxies: evolution; stars: formation; infrared; ultraviolet. 



1. Introduction 

Exploring the star formation history of galaxies is one of 
the most important topics in modem observational cosmol- 
ogy. Especially, the "true" absolute value of the cosmic star 
formation rate (hereafter SFR) has been of a central impor- 
tance to understand the formation and evolution of galaxies. 

However, it has been a difficult task for a long time be- 
cause of dust extinction. Active star formation (SF) is always 
accompanied by dust production through various dust grain 
formation processes related to the final stage of stellar evo- 
lution (e.g., Dwek and Scalo(1980); Dwek 1998; Nozawa 
et al. 2003; Takeuchi et al. 2005c; Asano et al. 2012). Ob- 
servationally, SFR of galaxies are, in principle, measured by 
the ultraviolet (UV) luminosity from massive stars because 
of their short lifetime (^ 10^ yr) compared with the age of 
galaxies or the Universe. However, since the UV photons are 
easily scattered and absorbed by dust grains, SFR of galaxies 
is always inevitably affected by dust produced by their own 
SF activity. Since the absorbed energy is re-emitted at wave- 
lengths of far-infrared (FIR), we need observations both at 
UV and FIR to have an unbiased view of their SF (e.g., Buat 
et al. 2005; Seibert et al. 2005; Cortese et al. 2006; Takeuchi 
et al. 2005a; Takeuchi et al. 2010; Haines et al. Bothwell et 
al. 2011;Haoetal. 2011). 
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After great effort of many researchers, the cosmic history 
of the SFR density is gradually converging at < z < 1. 
This "latter half" of the cosmic SFR is characterized by the 
rapid decline of the total SFR, especially the decrease of 
the contribution of dusty IR galaxies toward z = 0: While 
at z = 0, the contribution of the SFR hidden by dust is 
50-60 %, it becomes > 90 % at z = 1 (Takeuchi et al. 
2005a). This difference of decrease in SFR obtained from 
FUV and FIR has been already recognized in IR studies (e.g., 
Takeuchi et al. 2001a,b). Later works confirmed this gdusty 
era of the Universeh, and revealed that the dominance of the 
hidden SF continues even toward higher redshifts (z ~ 3) 
(e.g.. Murphy et al. 201 1; Cucciati et al. 2012). 

Then, a natural question arises: what does the different 
evolution at different wavelengths represent? To address this 
problem, it is very important to understand how we select 
sample galaxies and what we see in them. Each time we find 
some relation between different properties, we must under- 
stand clearly which is real (physical) and which is simply 
due to a selection effect. Some are detected at both bands, 
but some are detected only at one of the observed wave- 
length and appear as upper limits at the other wavelength. 
In previous studies it was often found that various claims 
were inconsistent with each other, mainly because they did 
not construct a well-controled sample of FUV and FIR se- 
lected galaxies. Recently, thanks to new large surveys, some 
attempts to explore the SFR distribution of galaxies in bivari- 
ate way have been made, through the "SFR function" (e.g.. 
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Buat et al. 2007, 2009, Haines et al. 2011; Bothwell et al. 
201 1). These works are based on the FUV and FIR LPs and 
their sum, but have not yet address their dependence on each 
other. To explore the bivariate properties of the SF in galax- 
ies further, the proper tool is the UV-IR bivariate luminosity 
function (BLF). 

However, constructing a BLF from two-wavelength data 
is not a trivial task. When we have a complete flux-limitec|3 
multiwavelength dataset, we can estimate a univariate lumi- 
nosity function (LF) at each band, but what we want to know 
is the dependence structure between luminosities at different 
bands. Mathematically, this problem is rephrased as follows: 
can we (re)construct a multivariate probability density func- 
tion (PDF) satisfying prescribed marginals? Although there 
is an infinite number of degrees of freedom to choose the 
original PDF, if we can model the dependence between vari- 
ables, we can construct such a bivariate PDF. A statistical 
tool for this problem is the so-called "copula" (see Sec. |2]for 
the definition). 

A copula has been extensively used in financial engineer- 
ing, for instance, but until recently there were very few ap- 
plications to astrophysical problems (e.g., Koen 2009; Ben- 
abed et al. 2009; ScheiTer et al. 2010; Sato et al. 2010, 
2011). Takeuchi (2010) introduced the copula to the esti- 
mation problem of a BLF. In this work, we apply the copula- 
based BLF estimation to the UV-IR two-wavelength dataset 
from z = to z = 1, using data from IRAS, AKAR^ 
Spitze^ and GALE:>^ 

The layout of this paper is as follows. We define a copula, 
especially the Gaussian copula, and formulate the copula- 
based BLF in Sec. |2] We then describe our UV and IR 
data in Sec. |3] In Sec. |4] we first formulate the likelihood 
function for the BLF estimation. Then we show the results, 
and discuss the possible interpretation of the evolution of the 
UV-IR BLF. SectionOis devoted to summary. 

Throughout the paper we will assume JImo = 0.3, f^Ao — 
0.7 and = 70 km s^^ Mpc^^. The luminosities are 
defined as i/L^ and expressed in solar units assuming Lq — 
3.83 X 10^3 ergs-^ 

2. The Bivariate Luminosity Function Based on the 

Copula 
2.1 Copula: general definition 

First, we briefly introduce a copula. A copula C(mi,M2) 
is defined as follows: 



we usually know the marginal DFs (or equivalently, PDFs) 
from the data. Then, the problem reduces to a statistical es- 
timation of a set of parameters to determine the shape of a 
copula C{ui ,U2)- In the form of the PDF, 



9{xi,x2) ^ dxidx2 fn^ifhi^^) 

= c[Fiixi),F2{x2)] fl{xi)f2{x2), 



(2) 



where fi{xi) and f2{x2) are the PDFs of _Fi(a;i) and 
F2{x2), respectively. On the second line, c(mi,M2) is re- 
ferred to as a differential copula. 
2.2 Gaussian copula 

Since the choice of copula is literally unlimited, we have 
to introduce a guidance principle. In many data analyses in 
physics, the most familiar measure of dependence might be 
the linear correlation coefficient p. Mathematically speak- 
ing, p depends not only on the dependence of two variables 
but also the marginal distributions, which is not an ideal 
property as a dependence measure. Even so, a copula hav- 
ing an explicit dependence on p would be convenient. In 
this work, we use a copula with this property, the Gaussian 
copula. 

One of the natural candidate with p may be a copula re- 
lated to a bivariate Gaussian DF (for other possibilities, see 
Takeuchi 2010). The Gaussian copula has an explicit depen- 
dence on a linear correlation coefficient by its construction. 
Let 



*i / ■^{x')dx\ 

<J —OQ 

1 



1p2ixi,X2; p) 



^(27r)2 det S 



exp 



(3) 



(4) 



(5) 



and 



4'2(a:i,a;2;p) 



ipi{x[, X2)dx[dx'2, (6) 



GixuX2)^C[Fiixi),F2{x2)] 



(1) 



where Fi{xi) and F2 {x2 ) are two univariate marginal cumu- 
lative distribution functions (DFs) and G{xi,X2) is a bivari- 
ate DF. We note that all bivariate DFs have this form and we 
can safely apply this method to any kind of bivariate DF es- 
timation problem (Takeuchi 2010). In various apphcations. 



where x = {xi, X2)'^ , S is a covariance matrix 



(7) 



^ Or any other observational condition. 

^ URL: |http://www.ir.isas.ac.jp/ASTRO-F/iiidex-e.html | 
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and superscript T stands for the transpose of a matrix or 
vector. 

Then, we define a Gaussian copula C'^ {ui,U2 \ p) as 



cG(^ii, p) = *2 [*r'(«i), *r'(«2); p] 



(8) 
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The density of C*^, c'^, is obtained as 

c°(ui,-u2;p) 

duidu2 

duidu2 

1 r 1 

Vd^^'^^X 2 



(9) 



where = *~^(u2)] and I stands for the 

identity matrix. 

2.3 Construction of the UV-IR BLF 

In this work, we define the luminosity at a certain wave- 
length band by L = vL^ {v is the corresponding frequency). 
Then the LF is defined as the number density of galax- 
ies whose luminosity lies between a logarithmic interval 

[logL, logL + dlog L\: 



dn 



dlogi 



(10) 



where we denote log a; = logj^g ^ ™d Inx = log^ For 
mathematical simplicity, we define the LF as being normal- 
ized, i.e.. 



j (f)^^\L)(l\ogL = 1 



(11) 



Hence, this corresponds to a PDF. We also define a cumula- 
tive LF as 



J log Lniin 



(12) 



where Lmin is the minimum luminosity of galaxies consid- 
ered. This corresponds to the DF. If we denote the univariate 
LFs as (t)i\Li) and ^^^(^2), then the BLF cjj'^^^Li, L2) is 
described by a differential copula c(ui , U2) as 



0(2)(£,,L2) = c[0W(Li)>«(L2) 
For the Gaussian copula, the BLF is obtained as 

0(2)(Li,L2;p) 



(13) 



1 fir 

Vd^'"'''^! 2 

x<t>['\L,)4'\L2) , 



I * 



where 



(14) 



(15) 



For the UV LF, we adopt the Schechter function 
(Schechter 1976). 



4')(L) = (lnlO) 0,1 (-^ 



L 



exp 



L 



(16) 



For the IR, we use the analytic form for the LF proposed by 
Saunders et al. (1990), 



^^2\L) 



L 



exp • 



2a2 



log 1 



L 



(17) 



We use the re-normalized version of eqs. (fTTt and (fT&t so 
that they can be regarded as PDFs, as mentioned above. 
2.4 Selection effects: another benefit of a copula BLF 

Another benefit of copula is that it is easy to incorpo- 
rate observational selection effects which always exist in any 
kind of astronomical data. In a bi(multi)variate analysis, 
there are two categories of observational selection effects. 

1) Truncation 

We do not know if a source would exist below a detec- 
tion limit. 

2) Censoring 

We know there is a source, but we have only an upper 
(sometimes lower) limit for a certain observable. 

As we mentioned above, we have to deal with both of these 
selection effects carefully to construct a BLF from data at 
the same time. It is terribly difficult to incorporate these 
effects by heuristic methods, particularly for a nonparametric 
methods (Takeuchi 2012, in preparation). In contrast, since 
we have an analytic form for the BLF, the treatment of upper 
limits is much more straightforward (Takeuchi 2010). We 
show how it is treated in the likelihood function in Sec. |4] 

3. Data 

3.1 UV-IR data construction 

We have constructed a dataset of galaxies selected at FUV 
and FIR by GALEX and IRAS for z = 0. At higher red- 
shifts, GALEX and Spitzer data are used for z = 0.7 and 
EIS (ESO Imaging Survey) and Spitzer for z = LO in the 
Chandra Deep Field South (CDFS). We explain the details 
of each sample in the following. 

In the Local Universe, we used the sample compiled by 
Buat et al. (2007). This sample was constructed based on 
IRAS all-sky survey and GALEX All-Sky Imaging Survey 
(AIS). This dataset consists of UV- and IR-selected sam- 
ples. The UV-selection was made by the GALEX FUV 
(A = 1530 A) band with FUV < 17 mag (hereafter, all 
magnitudes are presented in AB mag). This corresponds 
to the luminosity lower limit of ipuv > 10^ Lq. Red- 
shift information was taken from HyperLEDA (Paturel et 
al. 2003) and NED. The IR-selection is based on the IRAS 
PSCz (Saunders et al. 2000). The detection limit of the PSCz 
is 5*60 > 0.6 Jy. Redshifts of all PSCz galaxies were mea- 
sured, and most of their redshifts are z < 0.05. All UV 
fluxes were remeasured with the package we have developed 
(Iglesias-Paramo et al. 2006) to avoid the shredding of galax- 
ies. Details of the sample construction are explained in Buat 
et al. (2007). The number of the UV- and IR-selected sam- 
ples are 606 and 644, respectively. 

We also constructed a new, much larger sample of IR- 
selected galaxies by AKARIFIS All-Sky Survey. We started 



4 



T. T. TAKEUCHI et al: UV-IR BIVARIATE LUMINOSITY FUNCTION OF GALAXIES 



from the AKARI FIS bright source catalog (BSC) v. 1 from 
the AKARI all sky survey (Yamamura et al. 2010). This 
sample is selected at WIDE-S band (A — 90 /im) of the 
AKARI FIS (Kawada et al. 2007). The detection limit is 
5*90 > 0.2 Jy. We first selected AKARI sources in the 
SDSS footprints to have the same solid angle with a forth- 
coming corresponding UV-selected sample we are prepar- 
ing with GALEX-SDSSH Then, to have a secure sample 
of galaxies with redshift data, we made a cross match of 
AKAi?I sources with the Imperial IRAS-FSC Redshift Cata- 
logue (IIFSCz), a redshift catalog recently published (Wang 
and Rowan-Robinson 2009), with a search radius of 36 arc- 
sec. Since about 90 % of galaxies in the IIFSCz have spectro- 
scopic or photometric redshifts at Sqq > 0.36 Jy, the depth of 
the sample is defined by this matching. This determines the 
redshift range of this dataset, to be approximately z < 0.1. 
We measured the FUV and NUV flux densities with the same 
procedure as the IRAS-GALEX sample. The detection limits 
at FUV and NUV of this sample are 19.9 mag and 20.8 mag 
(Morrissey et al. 2007). A corresponding UV-selected sam- 
ple is under construction, hence we have only IR-selected 
sample. The number of galaxies is 3,891. For more informa- 
tion on this sample and properties of galaxies, see Sakurai et 
al. (2012). 

At higher-z, our samples are selected in the CDFS. 
GALEX observed this field at FUV and NUV (2300 A) as a 
part of its deep imaging survey. We restricted the field to the 
subfield observed by SpitzerMlPS as a part of the GOODS 
key program (Elbaz et al. 2007) to have the corresponding IR 
data. The area of the region is 0.068 deg^. Precise descrip- 
tion of our high-z samples are presented in Buat et al. (2009) 
and Burgarella et al. (2006). 

At z = 0.7, the NUV-band corresponds to the rest- 
frame FUV of the sample galaxies. We thus constructed 
the sample based on NUV-selection. Redshifts were taken 
from the COMBO- 17 survey (Wolf et al. 2004), and we 
have defined the z = 0.7 sample as those with redshifts 
of 0.6 < z < 0.8. Data reduction and redshift association 
are explained in Burgarella et al. (2006). We truncated the 
sample at NUV = 25.3 mag so that more than 90 % of the 
GALEX sources are identified in COMBO- 17 with R < 
24 mag. We set the MIPS 24 ^m upper limit as 0.025 mJy for 
undetected sources. For the IR-selected sample, we based on 
the GOODS SpitzerMlPS 24 fj,m sample and matched the 
GALEX and COMBO-17 sources. The sizes of the UV- and 
IR-selected samples are 340 and 470, respectively. 

In contrast to z = 0.7, since NUV of GALEX corresponds 
to 1 155 A in the rest frame of galaxies at z ~ 1.0, we cannot 
use NUV as the primary selection band as rest-frame FUV. 
Instead, we selected galaxies in the C/-band from the EIS 
survey (Arnouts et al. 2001) that covers the CDFS/GOODS 
field. We then cross-matched the U -band sources with the 
COMBO-17 sample to obtain redshifts. For the z = 1 
sample, the range of redshifts is 0.8 < z < 1.2. We set 
the limit atU = 24.3 mag to avoid source confusion. The IR 



^ This step is not a mandatory in this study, but we are planning to make an 
extension of this analysis with the UV-selected data, and this step will make 
the analysis easier with respect to the treatment of survey areas when the 
UV-selected data will be ready. 



flux densities were taken from SpitzerMlPS 24 fim data, and 
the same upper limit as the z = 0.7 sample is put for non- 
detections. Again the IR-selected sample was constructed 
from the GOODS and matched the GALEX and COMBO- 
17 sources, as the z ~ 0.7 sample. The sizes of the z = 1.0 
UV and IR-selected samples are 319 and 1033, respectively. 

The characteristics of these samples are summarized in 
Table [U 

3.2 Far-UV and total IR luminosities 

We are interested in the SF activity of galaxies and its 
evolution. Hence, luminosities of galaxies representative 
of SF activity would be ideal. As for the directly visible 
SF, obviously the UV emission is appropriate for this pur- 
pose. We define the FUV luminosity of galaxies Lpuv as 
Lpuv = J^Ly@FUV. For z = galaxies, FUV corresponds 
to 1530 A. At higher redshifts, as we explained, Lpuv is cal- 
culated from the NUV flux density at z — 0.7 and ?7-band 
flux density at z — 1.0, respectively. 

In contrast, at IR, the luminosity related to the SF activ- 
ity is the one integrated over a whole range of IR wave- 
lengths (A = 8-1000 ^m), Ltir, where the subscript 
TIR stands for the total IR. For the Local IRAS-GALEX 
sample, it is quite straightforward to define Ltir since the 
IRAS galaxies are selected at 60 fj,m. We adopted a for- 
mula Ltir — 2.5j^L^@60 fim. This rough approximation 
is, in fact, justified by spectral energy distributions (SEDs) 
of galaxies with ISO 160 /im observations (Takeuchi et al. 
2006). 

Also, for the AKARI-GALEX sample, we have excellent 
flux density data at FIR. We used the following TIR estima- 
tion formula from AKARI two wide bands (Takeuchi et al. 
2005b, 2010), 

Ltir Ai/(W/Z)£-S)L^(90 fim) 

+Aiy{WIDE-L)L^{UO urn), (18) 

where 

Aiy(WIDE-S) = 1.47 X 10^^ [Hzl 

(19) 

Av{WIDE-L) = 0.831 X 10^^ [Hz]. 

However, since deep FIR data of higher-z galaxies are not 
easily available to date, we should rely on the data taken by 
SpitzerMlPS 24 fjm. we use conversion formulae from MIR 
luminosity to Ltir: 

log Ltir [Lq] = 1.23 + 0.972 log Li5[Lq] , (20) 
logLTiR[LQ] = 2.27 + 0.707 log Li2[L0] 

+0.0140 (log Ll2[Lo])^ (21) 

which are an updated version of the formulae proposed by 
Takeuchi et al. (2005b) and also used by Buat et al. (2009). 
Here L15 and L12 are luminosities vL^@lbnm. and 12/im, 
i.e., 24/(1 + z) at z = 0.7 and 1.0, respectively. The es- 
timated Ltir slightly depends on which kind of conversion 
formula is used, but for our current purpose, it does not affect 
our conclusion and we do not discuss extensively here. In- 
tercomparison of the MIR-TIR conversion formulae can be 
found in Buat et al. (2009). 
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Table 1 


Sample description 




Redshift 


[range] 


Primary selection 


Number 


z = 


[0, 0.05] 


UV (GALEXFVV) 


606 






IR {IRAS 60 fim) 


644 




[0,0.1] 


IR (AKARI WIDE-S) 


3891 


z = 0.7 


[0.6,0.8] 


UV (GALEXNUY) 


340 






IR (Spitzer MIPS 24 /im) 


470 


z = 1.0 


[0.8, 1.2] 


UV (EIS u-band) 


319 






IR (Spitzer MIPS 24 /im) 


1033 



This situation will be greatly improved with Herschejl 
data. We will leave this as our future work with Herschel H- 
ATLAS (Bales et al. 2010) and GAMA (Galaxy And Mass 
Assembly: Driver et al. 201 1) data (Takeuchi et al. 2012, in 
preparation). 

4. Results and Discussion 
4.1 FUV and TIR univariate LFs 

In order to estimate the UV-IR BLF, first we have to ex- 
amine our setting for the FUV and TIR univariate LFs. The 
IRAS-GALEX sample, the validity of the Local univariate 
LFs are already proved [see Fig. 3 of Buat et al. (2007)]. 
Hence, we can safely use the analytic formulae of FUV and 
TIR LFs at z = [eqs. ^ and 

We use the Schechter parameters presented by Wyder et 
al. (2005) for GAL£X FUV: (ai,L*i,0,i) = (1.21,1.81 x 
109/1-2 [Lq],1.35 X IQ-^h^ [Mpc'3]). For the TIR, 
we used the parameters (a2, i*2, </'*2, f) = (1.23,4.34 x 
lO^h-^ [L0],2.34x IQ-^h^ []VIpc-^j, 0.724) (Takeuchi et 
al. 2003) obtained from the IRAS PSCz galaxies (Saunders 
et al. 2000), and multiplied with 2.5 to convert the 60- 
^m LF to the TIR LF 

For higher redshifts, ideally we should estimate the uni- 
variate FUV and TIR LFs simultaneously with the BLF esti- 
mation. It is, however, quite difficult for our cuiTent samples 
because of the limited number of galaxies. We instead used 
the LF parameters at z — 0.7 and 1.0 obtained by previ- 
ous studies on univariate LFs and modeled the FUV and TIR 
univariate LFs at these redshifts and examined their validity 
with nonparametric LFs estimated from the data. We use the 
parameters compiled by Takeuchi et al. (2005a). Parameters 
of the evolution of the TIR LF are obtained by approximating 
the evolution in the form 



f(^) 



(22) 



where (f'2^l{z) is the local functional form of the TIR LF. 
Le Floc'h et al. (2005) assumed a power-law form for the 
evolution functions as 



f(z) = (l + z)«, g{z) = {l + z) 



(23) 



and obtained *p = 0.7 and = 3.2, with a remaining 
constant. The Schechter parameters at z = 0.7 and 1.0 for 



the FUV LF are directly estimated by Arnouts et al. (2005) 
and we adopt their values (see Table 1 of Takeuchi et al. 
2005a). 

Then we estimated the FUV and TIR LFs with the step- 
wise maximum likelihood method and the variant of C~ 
method from our CDFS multiwavelength data (for the es- 
timation method, see, e.g., Takeuchi et al. 2000; Johnston 
2011, and references therein). The obtained univariate LFs 
at z — 0.7 and 1.0 are presented in Figs. [T] and |2] We also 
show the analytic model LFs in these figures. We note that a 
well-known large density enhancement locates in the CDFS 
at z ~ 0.7 (e.g., Salimbeni et al. 2009), and thus we have 
renormalized the LFs to remove the effect of the overdensity. 
In the following analysis, we normalize the univariate LFs 
according to eq. ( fTTT i so that we can treat the univariate LFs 
as PDFs, hence this does not affect the following analysis at 
all. 

Because of the small sample size, the LF shape does not 
perfectly agree with the supposed functional forms, but the 
nonparametric LFs are acceptably similar to the analytic 
functions both for the FUV and TIR at each redshift. We 
stress that the analytic functions are not the fit to the data, 
but estimated from other studies. This implies that the es- 
timated evolutionary parameters of the LFs work generally 
well. Thus, we can use the higher-z univariate LFs as the 
marginal PDFs for the estimation of the UV-IR BLFs. 
4.2 Copula likelihood for the BLF estimation 

By using the estimated univariate FUV and TIR LFs 
as given marginals, we can estimate only one parameter, 
the linear correlation p by maximizing the likelihood func- 
tion. The structure of a set of two-band selected data is 
(^fuV'Juv'-^tiR'Jir)'**; = l...nk. Here, jband (band: 
UV or IR) stands for the upper limit flag: jband — 0: detec- 
tion and jband — upper limit. Another index k is the 
indicator of the selected band, i.e., k = 1 means a sample 
galaxy is selected at UV and fc = — 1 means it is selected at 



' URL: |http://hersctiel.esac.esa.int/ 1 
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IR. The likelihood function is as follows. 



ln£ (Lp'^uy,i^fjj^| ife 1, • • - rife, fc = 1, -l) 
"fc ( 



where p'*'^' (ipyy, -^tir) '■^^ probability for ikth galaxy 
to be detected at both bands and to have luminosities Lpuy 
and L^jj^, 

„det f fik T ik 



detfrik jik \ 



E E^-[/^'(^PUV,i^TlH)]^^"^"^^^^^"^'^^^ 

, _ r 1 : UV scl Zfc^l L 
'^"l -1 : IR acl 



(2) (TL 



ln[p"^^"M^FW,iT\R)] 
-ltt[p"L^'^(i-uy,L-j,,)] 



(^FUVi -^TIR; P) d-C'TIRdi^FUV 

(25) 



„UL:UV { rik 



(l-fc)(-3 " ) 



(Lpuy, is the probability for ifcth galaxy to 

(24)^^ detected at IR band and have a luminosity ^xiii' ^'^^ 
detected at UV band and only an upper limit ipuv jk 
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Fig. 3. Schematic BLF. (1) Diagonal: The energy from SF is emitted 
equally at UV and IR with any SF activity. (2) Upward: The more 
active the SF in a galaxy is, the more luminous at the IR (dusty SF). 
(3) Downward: The more active the SF is, the more luminous at the UV 
(gtransparenth SF). 
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(26) 

Fuv ^tir) probability for ikth galaxy 

to be detected at UV band and have a luminosity L^-^jy, but 
not detected at IR band and only an upper limit i^iR jk 
available, 

UL: IR (r ik T Ik \ 
f l-^FUV'^TIRj 



TIR 



/ / 0^^^ (iFUV>^TIR;P)diTIRdiFUV 

•^i™v(^.J-^0 

(27) 

The denominator in eq. dZSl l is introduced to take into ac- 
count the truncation in the data by observational flux selec- 
tion limits at both bands (e.g., Sandage et al. 1979; Johnston 
2011). We should also note that it often happens that the 
same galaxies are included in both the UV- and IR-selected 
sample. In such a case, we should count the same galaxies 
only once to avoid double counting of them. Practically, such 
galaxies are included in any of the samples, because they are 
detected at both bands and are symmetric between UV- and 
IR-selections. 

In the estimation procedure, we estimated the univariate 
LFs and their evolutions first, and then used these parameters 



to estimate the BLF and its evolution. One might wonder if 
this compromising method would introduce some bias in the 
estimation. Since we fix the form of the maximum likelihood 
estimator eq. (l24l i. this two-step estimation instead of simul- 
taneous estimation does not bias the result, especially the de- 
pendence structure between UV and IR, unless the assumed 
univariate LP shape would be significantly different from the 
real one. We have already seen that the stepwise maximum 
likelihood estimation gave nonparametric LFs which reason- 
ably agree with the assumed Schechter or Saunders function. 
Then, we can safely rely on the result for further discussion. 

However, we should note that the two-step estimation 
gives conditional errors of the parameters only for each step, 
not the marginal ones. Then, the errors for each parame- 
ters are significantly underestimated. If we want to discuss 
the error values more precisely, we need a larger dataset and 
must use the simultaneous estimation method. This will be 
done in the future works. 
4.3 The BLF and its evolution 

Using the Gaussian copula, now we can estimate the BLF. 
The visible and hidden SFRs should be directly reflected 
to this function. Dust is produced by SF activity, but also 
destroyed by SN blast waves as a result of the SF. Many 
physical processes are related to the evolution of the dust 
amount. Thus, first of all, we should describe statistically 
how it evolved, as stated in Introduction. 

First, we summarize how to interpret the UV-IR BLF 
schematically in Fig. [3] First, we see the case that the ridge 
of the BLF is straight and diagonal [see (1) in Fig.|3]. This 
means that the energy from SF is emitted equally at UV and 
IR with any SF activity. If the relation is diagonal but has an 
offset horizontally or vertically, this suggests that a constant 
fraction of energy is absorbed by dust and reemitted. Second, 
if the ridge is curved upward, it means that the more active 
the SF in a galaxy is, the more luminous at the IR [dusty SF: 

(2) in Fig.O. Third, if the ridge is bent downward, the more 
active the SF is, the more luminous at UV ['transparent' SF: 

(3) in Fig. El . 

Now we show the estimated BLFs in Figs. |4]-[7] In the 
Local Universe, the BLF is quite well constrained. The esti- 
mated correlation coefficient p is very high: p = 0.95 ± 0.04 
for IRAS-GALEX and p = 0.95 ± 0.006 for AKARI- 
GALEX datasets. The apparent scatter of the ipuv-iTiR 
is found to be due to the nonlinear shape of the ridge of 
the BLF. This bent shape of the BLF was implied by pre- 
ceding studies (Martin et al. 2005), and we could quantify 
this feature. The copula BLF naturally reproduced it. For 
the AKARI-GALEX sample, only the IR-selected sample is 
available at this moment. This, however, does not bias the 
BLF estimation since the information from censored data 
points are incorporated properly in eq. (1241 . This can be di- 
rectly tested, for example, by using only one of the UV and 
IR-selected data for the z — BLF estimation with IRAS- 
GALEX dataset. Both one-band estimations with UV and 
IR-selected data yield p = 0.95 ± 0.07. Obviously the error 
becomes larger because of fewer data, but the estimate itself 
remain unchanged. 

At higher redshifts (z — 0.7-1.0), the linear correlation re- 
mains tight (p ~ 0.91 ± 0.06 at z = 0.7 and p ~ 0.86 ± 0.05 
at z = 1.0), even though it is difficult to constrain the 
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Fig. 4. The BLF of galaxies from IRAS and GALEX at z = 0. Contours 
are the analytic model constructed by a Gaussian copula and univariate 
FUV and TIR LFs. Open squares represent the UV-selected sample 
from GALEX, while open triangles are the IR-selected ones from IRAS. 
Squares and triangles with arrows mean that they are upper limits at FUV 
and FIR, respectively. 
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FUV lumi nosity ^puv L^oJ 

Fig. 5. The same as Fig. |4] but from AKARI and GALEX. Open trian- 
gles represent the IR-selected samples from AKARI. Since there is no 
UV-selected sample in this figure, we only show the IR-selected ones. 



low-luminosity end from the data in this analysis (Spitzer- 
GALEX in the CDFS). The distribution of upper limits in 
Figs. |6] and Q looks different from that in Fig. |4] Since we 
have restricted the redshifts of galaxies in these samples, the 
redshift restriction gives approximately constant luminosity 
limits both at FUV and FIR. This gives the "L-shaped" dis- 
tribution of upper limits seen in these figures. 

Though the whole shape cannot be perfectly determined 
by the current data, we find that p in the copula LF is high 
and remarkably stable with redshifts (from 0.95 at z = to 



Fig. 6. The same as Fig. |4] but from Spitzer and GALEX at z = 0.7. 
Symbols are essentially the same as in Fig.|4] but IR-selected samples are 
from SpitzeiMlPS. 



10" 




FUV luminosity L^^y [Lq] 



Fig. 7. The same as Fig. |4]but from Spitzer and GALEX at 2 = 1.0. Again 
symbols are essentially the same as in Fig.|4] but IR-selected samples are 
from Spitzei/MIPS. 

0.85 at z = 1.0). This impHes the evolution of the UV-IR 
bivariate LF is mainly due to the different evolution of the 
univariate LFs, and may not be controled by the evolution of 
the dependence structure. 

Forthcoming better data in the future will improve various 
aspects of the BLF estimation. First, Herschel, ALMA, and 
will provide us with direct estimations of Ltir of 
galaxies especially at high redshifts, since they will detect 
these galaxies at longer wavelengths than 24 /im. This allows 
us to avoid a possible bias of the estimated Ltir caused by 



URL: http://www.iT.isas.jaxa.jp/SPICA/SPICA_HP/index_EngIish.htmll 
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the 24 /xm-based extrapolation of the SED. 

Second, when we have deeper data at UV and IR, we 
will be able to extend the luminosity range toward much 
lower luminosities both at FUV and FIR. Then, the faint-end 
structure of the BLF will be much more tightly determined. 
Further, it is not only for the determination of the faint end, 
but deeper data will give some insights for the choice of 
the copula functional form. In this study, we did not try to 
examine if the Gaussian copula would be a proper choice as 
a model of the UV-IR BLF. At z = 0, when UV-selected 
data comparably large as the AKARI sample are ready, we 
will be able to specify (a family of) copulas appropriate for 
this analysis. At higher redshifts, our current data are not 
deep enough to have any claim on the faint-end dependence 
structure of the BLF, but again Herschel, ALMA, and SPICA 
data together with ground-based optical or JWST ones will 
enables us to examine which copula would be appropriate 
to describe the BLF, and constrain the evolution of the BLF 
along the whole history of galaxy evolution in the Universe. 

5. Conclusion 

To understand the visible and hidden SF history in the 
Universe, it is crucial to analyze multiwavelength data in a 
unified manner. The copula is an ideal tool to combine two 

marginal univariate LFs to construct a bivariate LFs. It is 
straightforward to extend this method to multivariate DFs. 

1) The Gaussian copula LF is sensitive to the linear corre- 
lation parameter p. 

2) Even so, p in the copula LF is remarkably stable with 
redshifts (from 0.95 at z = to 0.85 at z = 1.0). 

3) This implies the evolution of the UV-IR BLF is mainly 
due to the different evolution of the univariate LFs, and 
may not be controled by the dependence structure. 

4) The nonlinear structure of the BLF is naturally repro- 
duced by the Gaussian copula. 

The data from Herschel, ALMA, and SPICA data will im- 
prove the estimates drastically, and we expect to specify the 
full evolution of the UV-IR BLF in the Universe. We stress 
that the copula will be a useful tool for any other kind of bi- 
(muM-) variate statistical analysis. 
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